#1
th = seq(0,1,length=100)

plot(th,dbeta(th,2,2),type="l",xlab=expression(pi),ylab="",ylim=c(0.0,6.5),col="red")
par(new=T)
plot(th,dbeta(th,2+22,2+35),type="l",xlab="",ylab="",col="blue",axes=F,ylim=c(0.0,6.5))

plot(th,dbeta(th,5,10),type="l",xlab=expression(pi),ylab="",ylim=c(0.0,7),col="red")
par(new=T)
plot(th,dbeta(th,5+22,10+35),type="l",xlab="",ylab="",col="blue",axes=F,ylim=c(0.0,7))

plot(th,dbeta(th,100,3),type="l",xlab=expression(pi),ylab="",ylim=c(0.0,30),col="red")
par(new=T)
plot(th,dbeta(th,100+22,3+35),type="l",xlab="",ylab="",col="blue",axes=F,ylim=c(0.0,30))

#2
hpd(qbeta, shape1=24, shape2=37,conf=0.95)
hpd(qbeta, shape1=27, shape2=45,conf=0.95)
hpd(qbeta, shape1=122, shape2=38,conf=0.95)
24/(24+37)
27/(27+45)
122/(122+38)

#3 want to get f(mu|Y)
x = rnorm(20,10,4)
Xt <- rchisq(1000,19)

sig1 = 4
mut1 = rnorm(1000,mean(x),sqrt(sig1/20)) # known sig, f(mu) proportional to 1, get f(mu|Y)
emp.hpd(mut1,conf=0.95)
quantile(mut1,c(0.025,0.975))

sig2 = 19*var(x)/Xt
mut2 = rnorm(1000,mean(x),sqrt(sig2/20)) # this isn't correct, get f(mu|sig^2,Y)
emp.hpd(mut2,conf=0.95)

mut3 = rt(1000,19)*sqrt(var(x)/20)+mean(x) # sig unknown, f(mu,sig^2) proportional to 1/sig^2, get f(mu|Y)
emp.hpd(mut3,conf=0.95)
quantile(mut2,c(0.025,0.975))

# TODO try quantile

#4
muA <- rt(1000,6)*sqrt(52.52/7)+80.63
muB <- rt(1000,6)*sqrt(243/7)+50.18
diffs <- muA - muB
hist(diffs,main="",ylab="",xlab=expression(mu[A]-mu[B]))
quantile(diffs,c(0.025,0.975))
